Chapter 4 Diversity analysis
genome_counts_filt <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
select(where(~ sum(.) > 0)) %>%
select(-AC85) %>%
rownames_to_column(., "genome")
genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt$genome)
table <- genome_counts_filt %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character,as.numeric) %>%
rownames_to_column(., "sample")
master_table <- sample_metadata %>%
mutate(sample=Tube_code) %>%
mutate(Tube_code=str_remove_all(Tube_code, "_a")) %>%
filter(Tube_code %in% table$sample) %>%
mutate(treatment = sub("^\\d+_", "", time_point))%>%
left_join(., table, by=join_by("Tube_code" =="sample"))
sample_metadata <- master_table %>%
select(1:13)
genome_counts_filt <- master_table %>%
select(12,14:323) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character,as.numeric) %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
dplyr::select(where(~ !all(. == 0)))%>%
rownames_to_column(., "genome")
genome_counts_filtering <- master_table %>%
select(12,14:323) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character,as.numeric) %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
dplyr::select(where(~ !all(. == 0)))
genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt$genome)
#match_data(genome_counts_filtering,genome_tree)
genome_metadata <- genome_metadata %>%
filter(genome %in% genome_counts_filt$genome)4.1 Calculate diversities
4.1.1 Alpha diversity
# Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 0) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(richness = 1) %>%
rownames_to_column(var = "sample")
neutral <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(neutral = 1) %>%
rownames_to_column(var = "sample")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, tree = genome_tree) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(phylogenetic = 1) %>%
rownames_to_column(var = "sample")
genome_gifts1 <- genome_gifts[genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]
genome_counts_filt1 <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts1),]
genome_counts_filt1 <- genome_counts_filt1 %>%
remove_rownames() %>%
column_to_rownames(var = "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
rownames_to_column(., "genome")
# Merge all metrics
alpha_div <- richness %>%
full_join(neutral, by = join_by(sample == sample)) %>%
full_join(phylogenetic, by = join_by(sample == sample))4.1.2 Beta diversity
beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 0)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1)
beta_q1p <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, tree = genome_tree)
genome_gifts1 <- genome_gifts[genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]
genome_counts_filt1 <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts1),]4.2 Does captivity change the microbial community?
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(time_point %in% c("Acclimation", "Wild") ) 4.2.1 Alpha diversity
label_map <- c(
"Cold_wet" = "Cold population",
"Hot_dry" = "Warm-population",
"richness" = "Species Richness",
"neutral" = "Neutral Diversity",
"phylogenetic" = "Phylogenetic Diversity")
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(
metric = factor(metric, levels = c("richness", "neutral", "phylogenetic")),
time_point = factor(time_point, levels = c("Wild", "Acclimation"))) %>%
filter(time_point %in% c("Wild", "Acclimation")) %>%
ggplot(aes(y = value, x = time_point, color=Population, fill=Population)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
facet_grid(metric ~ Population, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ Population*time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(9.4635) ( log )
Formula: richness ~ Population * time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
518.4 530.2 -253.2 506.4 47
Scaled residuals:
Min 1Q Median 3Q Max
-2.07994 -0.34567 0.00022 0.55317 1.22959
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0.1126 0.3355
Number of obs: 53, groups: individual, 27
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.7707 0.1208 31.217 < 2e-16 ***
PopulationHot_dry 0.6756 0.2006 3.368 0.000757 ***
time_pointWild 0.4312 0.1223 3.525 0.000424 ***
PopulationHot_dry:time_pointWild -0.5331 0.2048 -2.604 0.009221 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) PpltH_ tm_pnW
PpltnHt_dry -0.602
time_pntWld -0.557 0.335
PpltnHt_:_W 0.338 -0.526 -0.601$emmeans
time_point emmean SE df asymp.LCL asymp.UCL
Acclimation 4.11 0.100 Inf 3.91 4.31
Wild 4.27 0.099 Inf 4.08 4.47
Results are averaged over the levels of: Population
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Acclimation - Wild -0.165 0.102 Inf -1.615 0.1063
Results are averaged over the levels of: Population
Results are given on the log (not the response) scale. $emmeans
Population emmean SE df asymp.LCL asymp.UCL
Cold_wet 3.99 0.101 Inf 3.79 4.18
Hot_dry 4.40 0.138 Inf 4.12 4.67
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Cold_wet - Hot_dry -0.409 0.171 Inf -2.396 0.0166
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale. Neutral
Model_neutral <- nlme::lme(fixed = neutral ~ time_point*Population, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
402.7832 414.1342 -195.3916
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 5.51778 10.54924
Fixed effects: neutral ~ time_point * Population
Value Std.Error DF t-value p-value
(Intercept) 19.37729 2.883717 25 6.719554 0.0000
time_pointWild 17.27881 3.578682 24 4.828260 0.0001
PopulationHot_dry 24.71329 4.905493 25 5.037882 0.0000
time_pointWild:PopulationHot_dry -22.28917 6.126768 24 -3.637998 0.0013
Correlation:
(Intr) tm_pnW PpltH_
time_pointWild -0.642
PopulationHot_dry -0.588 0.377
time_pointWild:PopulationHot_dry 0.375 -0.584 -0.632
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.5974563 -0.5606149 0.1376821 0.5636829 1.6703831
Number of Observations: 53
Number of Groups: 27 $emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 31.7 2.45 25 26.7 36.8
Wild 37.9 2.43 24 32.9 42.9
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Wild -6.13 3.06 24 -2.002 0.0567
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment $emmeans
Population emmean SE df lower.CL upper.CL
Cold_wet 28.0 2.21 24 23.5 32.6
Hot_dry 41.6 3.09 24 35.2 48.0
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Cold_wet - Hot_dry -13.6 3.8 24 -3.568 0.0016
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment Phylogenetic
Model_phylo <- nlme::lme(fixed = phylogenetic ~ time_point*Population, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
193.7567 205.1077 -90.87837
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.7860836 1.191838
Fixed effects: phylogenetic ~ time_point * Population
Value Std.Error DF t-value p-value
(Intercept) 5.343318 0.3453897 25 15.470404 0.0000
time_pointWild -0.135493 0.4048212 24 -0.334698 0.7408
PopulationHot_dry 1.119702 0.5880336 25 1.904146 0.0685
time_pointWild:PopulationHot_dry -0.602016 0.6924897 24 -0.869350 0.3933
Correlation:
(Intr) tm_pnW PpltH_
time_pointWild -0.608
PopulationHot_dry -0.587 0.357
time_pointWild:PopulationHot_dry 0.355 -0.585 -0.596
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.2027374 -0.4577650 0.0610368 0.3504732 1.7675077
Number of Observations: 53
Number of Groups: 27 $emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 31.7 2.45 25 26.7 36.8
Wild 37.9 2.43 24 32.9 42.9
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Wild -6.13 3.06 24 -2.002 0.0567
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment $emmeans
Population emmean SE df lower.CL upper.CL
Cold_wet 28.0 2.21 24 23.5 32.6
Hot_dry 41.6 3.09 24 35.2 48.0
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Cold_wet - Hot_dry -13.6 3.8 24 -3.568 0.0016
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment 4.2.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(time_point %in% c("Acclimation", "Wild")) %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(time_point %in% c("Acclimation", "Wild"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.03445 0.034446 4.6041 999 0.036 *
Residuals 51 0.38156 0.007482
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(richness ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 0.5784094 | 0.03513500 | 2.236493 | 0.001 |
| Population | 1 | 2.8737192 | 0.17456169 | 11.111600 | 0.001 |
| time_point:Population | 1 | 0.3378127 | 0.02052015 | 1.306196 | 0.001 |
| Residual | 49 | 12.6725436 | 0.76978316 | NA | NA |
| Total | 52 | 16.4624849 | 1.00000000 | NA | NA |
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(richness),]
pairwise <- pairwise.adonis(richness, subset_meta_arrange$Population, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_wet vs Hot_dry | 1 | 2.87727 | 10.8015 | 0.1747774 | 0.001 | 0.001 | ** |
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Wild vs Acclimation | 1 | 0.5784094 | 1.857135 | 0.035135 | 0.023 | 0.023 | . |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.01130 0.0113018 1.2646 999 0.259
Residuals 51 0.45578 0.0089369 adonis2(neutral ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 0.8535852 | 0.05306573 | 3.688514 | 0.001 |
| Population | 1 | 3.3711165 | 0.20957575 | 14.567274 | 0.001 |
| time_point:Population | 1 | 0.5212922 | 0.03240772 | 2.252609 | 0.001 |
| Residual | 49 | 11.3394385 | 0.70495080 | NA | NA |
| Total | 52 | 16.0854325 | 1.00000000 | NA | NA |
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(neutral),]
pairwise <- pairwise.adonis(neutral, subset_meta_arrange$Population, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_wet vs Hot_dry | 1 | 3.369473 | 13.51397 | 0.2094736 | 0.001 | 0.001 | ** |
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Wild vs Acclimation | 1 | 0.8535852 | 2.858015 | 0.05306573 | 0.002 | 0.002 | * |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00207 0.0020675 0.125 999 0.735
Residuals 51 0.84374 0.0165440 adonis2(phylo ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 0.2282721 | 0.07008479 | 4.486412 | 0.001 |
| Population | 1 | 0.3478468 | 0.10679697 | 6.836508 | 0.001 |
| time_point:Population | 1 | 0.1878081 | 0.05766139 | 3.691140 | 0.003 |
| Residual | 49 | 2.4931580 | 0.76545685 | NA | NA |
| Total | 52 | 3.2570850 | 1.00000000 | NA | NA |
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(phylo),]
pairwise <- pairwise.adonis(phylo, subset_meta_arrange$Population, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_wet vs Hot_dry | 1 | 0.3462509 | 6.066576 | 0.106307 | 0.002 | 0.002 | * |
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Wild vs Acclimation | 1 | 0.2282721 | 3.84371 | 0.07008479 | 0.004 | 0.004 | * |
dbRDA
#Richness
cca_ord <- capscale(formula = richness ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta$time_pointWild" = "Wild",
"subset_meta$PopulationHot_dry"="Warm population",
"subset_meta$time_pointWild:subset_meta$PopulationHot_dry"="Interaction Wam population")
beta_richness_rda <-CAP_df %>%
group_by(Population, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta$time_pointWild" = "Wild",
"subset_meta$PopulationHot_dry"="Warm population",
"subset_meta$time_pointWild:subset_meta$PopulationHot_dry"="Interaction Wam population")
beta_neutral_rda <- CAP_df %>%
group_by(Population, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50"))+
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta$time_pointWild" = "Wild",
"subset_meta$PopulationHot_dry"="Warm population",
"subset_meta$time_pointWild:subset_meta$PopulationHot_dry"="Interaction Wam population")
beta_phylogenetic_rda<- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()ggarrange(beta_richness_rda, beta_neutral_rda, beta_phylogenetic_rda, ncol=3, nrow=1, common.legend = TRUE, legend="right")
4.2.3 Cold population
4.2.3.1 Structural zeros
wild_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point == "Wild") %>%
dplyr::select(sample) %>%
pull()
acclimation_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point == "Acclimation") %>%
dplyr::select(sample) %>% pull()
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",subset_meta$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_wild = all(c_across(all_of(wild_samples)) == 0)) %>%
mutate(all_zeros_acclimation= all(c_across(all_of(acclimation_samples)) == 0)) %>%
mutate(average_wild = mean(c_across(all_of(wild_samples)), na.rm = TRUE)) %>%
mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>%
filter(all_zeros_wild == TRUE || all_zeros_acclimation==TRUE) %>%
mutate(present = case_when(
all_zeros_wild & !all_zeros_acclimation ~ "acclimation",
!all_zeros_wild & all_zeros_acclimation ~ "wild",
!all_zeros_wild & !all_zeros_acclimation ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "acclimation", average_wild, average_acclimation)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Wild
structural_zeros %>%
filter(present=="wild") %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) # A tibble: 5 × 2
# Rowwise:
phylum sample_count
<chr> <int>
1 p__Bacillota_A 7
2 p__Bacillota 1
3 p__Bacillota_B 1
4 p__Bacteroidota 1
5 p__Spirochaetota 1Acclimation
# A tibble: 14 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_12:bin_000001 p__Pseudomonadota c__Gammaproteobacteria o__Pseudomonadales f__Moraxellaceae g__Acinetobacter s__Acinetobacter jo…
2 AH1_2nd_15:bin_000001 p__Pseudomonadota c__Alphaproteobacteria o__Rhizobiales f__Rhizobiaceae g__Agrobacterium s__Agrobacterium tu…
3 AH1_2nd_17:bin_000010 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Rikenellaceae g__Mucinivorans s__
4 AH1_2nd_17:bin_000038 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides_G s__
5 AH1_2nd_20:bin_000014 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Citrobacter s__Citrobacter port…
6 AH1_2nd_20:bin_000061 p__Bacillota c__Bacilli o__Lactobacillales f__Enterococcaceae g__Enterococcus s__
7 AH1_2nd_20:bin_000073 p__Bacillota c__Bacilli o__Lactobacillales f__Enterococcaceae g__Enterococcus s__Enterococcus sp0…
8 LI1_2nd_10:bin_000001 p__Bacillota c__Bacilli o__Lactobacillales f__Enterococcaceae g__Enterococcus_A s__Enterococcus_A r…
9 LI1_2nd_10:bin_000017 p__Chlamydiota c__Chlamydiia o__Chlamydiales f__Chlamydiaceae g__ s__
10 LI1_2nd_2:bin_000039 p__Bacillota_A c__Clostridia o__TANB77 f__CAG-508 g__ s__
11 LI1_2nd_7:bin_000045 p__Bacillota_A c__Clostridia o__Oscillospirales f__Acutalibacteraceae g__ s__
12 LI1_2nd_7:bin_000074 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Escherichia s__Escherichia coli
13 LI1_2nd_8:bin_000030 p__Actinomycetota c__Actinomycetia o__Mycobacteriales f__Mycobacteriaceae g__Corynebacterium s__
14 LI1_2nd_8:bin_000079 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Citrobacter_A s__Citrobacter_A am…Community elements differences in structural zeros
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>4.2.3.2 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point %in% c("Wild", "Acclimation") )%>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "time_point",
rand_formula = "(1|individual)",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output2$res %>%
dplyr::select(taxon, lfc_time_pointWild, p_time_pointWild) %>%
filter(p_time_pointWild < 0.05) %>%
dplyr::arrange(p_time_pointWild) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_time_pointWild)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_pointWild, y=forcats::fct_reorder(genome,lfc_time_pointWild), fill=phylum)) + #forcats::fct_rev()
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
# coord_flip()+
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))
Phyla of the significant MAGs
ancombc_rand_table_mag%>%
filter(lfc_time_pointWild<0) %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) phylum sample_count
1 Bacillota_A 5
2 Bacillota 4
3 Bacteroidota 1
4 Campylobacterota 1
5 Cyanobacteriota 1 phylum genus species
1 Bacillota_A
2 Bacillota Coprobacillus
3 Bacillota Mycoplasmoides
4 Cyanobacteriota Scatousia
5 Bacteroidota Bacteroides
6 Bacillota NSJ-61
7 Bacillota_A
8 Campylobacterota Helicobacter_J
9 Bacillota Clostridium_AQ
10 Bacillota_A CHH4-2
11 Bacillota_A Enterocloster
12 Bacillota_A UBA866 ancombc_rand_table_mag%>%
filter(lfc_time_pointWild<0) %>%
count(genus, name = "sample_count") %>%
arrange(desc(sample_count)) genus sample_count
1 2
2 Bacteroides 1
3 CHH4-2 1
4 Clostridium_AQ 1
5 Coprobacillus 1
6 Enterocloster 1
7 Helicobacter_J 1
8 Mycoplasmoides 1
9 NSJ-61 1
10 Scatousia 1
11 UBA866 14.2.3.3 Differences in the functional capacity
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Wild"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.2.3.3.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.335 0.0324
2 Wild 0.350 0.0237MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94347, p-value = 0.07144
Wilcoxon rank sum exact test
data: value by treatment
W = 114, p-value = 0.2066
alternative hypothesis: true location shift is not equal to 0element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-Wild)%>%
mutate(group_color = ifelse(Difference <0, "Wild","Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=c('#008080', '#00808050')) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")
4.2.3.3.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.330 0.0320
2 Wild 0.346 0.0196MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94507, p-value = 0.07998
Wilcoxon rank sum exact test
data: value by treatment
W = 100, p-value = 0.08313
alternative hypothesis: true location shift is not equal to 0function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))| Code_function | Acclimation | Wild | Function | Difference | treatment |
|---|---|---|---|---|---|
| B03 | 0.1097356 | 0.1256091 | Amino acid derivative biosynthesis | -0.0158735 | Wild |
| D03 | 0.2921075 | 0.3442827 | Sugar degradation | -0.0521752 | Wild |
4.2.4 Warm population
4.2.4.1 Structural zeros
wild_samples_warm <- sample_metadata %>%
filter(Population == "Hot_dry" & time_point == "Wild") %>%
dplyr::select(sample) %>%
pull()
acclimation_samples_warm <- sample_metadata %>%
filter(Population == "Hot_dry" & time_point == "Acclimation") %>%
dplyr::select(sample) %>% pull()
sample_sub <- sample_metadata %>%
filter(Population == "Hot_dry" & treatment %in% c("Acclimation", "Wild"))
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",sample_sub$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_wild_warm = all(c_across(all_of(wild_samples_warm)) == 0)) %>%
mutate(all_zeros_acclimation_warm= all(c_across(all_of(acclimation_samples_warm)) == 0)) %>%
mutate(average_wild_warm = mean(c_across(all_of(wild_samples_warm)), na.rm = TRUE)) %>%
mutate(average_acclimation_warm = mean(c_across(all_of(acclimation_samples_warm)), na.rm = TRUE)) %>%
filter(all_zeros_wild_warm == TRUE || all_zeros_acclimation_warm==TRUE) %>%
mutate(present = case_when(
all_zeros_wild_warm & !all_zeros_acclimation_warm ~ "acclimation",
!all_zeros_wild_warm & all_zeros_acclimation_warm ~ "wild",
!all_zeros_wild_warm & !all_zeros_acclimation_warm ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "acclimation", average_wild_warm, average_acclimation_warm)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Wild
structural_zeros %>%
filter(present=="wild") %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) # A tibble: 7 × 2
# Rowwise:
phylum sample_count
<chr> <int>
1 p__Bacillota_A 11
2 p__Bacillota 4
3 p__Pseudomonadota 4
4 p__Bacillota_C 2
5 p__Bacteroidota 1
6 p__Desulfobacterota 1
7 p__Fusobacteriota 1Acclimation
# A tibble: 32 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_10:bin_000010 p__Cyanobacteriota c__Vampirovibrionia o__Gastranaerophilales f__Gastranaerophilaceae g__CALURL01 s__
2 AH1_2nd_12:bin_000001 p__Pseudomonadota c__Gammaproteobacteria o__Pseudomonadales f__Moraxellaceae g__Acinetobacter s__Acinetoba…
3 AH1_2nd_14:bin_000015 p__Bacillota_A c__Clostridia o__Christensenellales f__UBA3700 g__ s__
4 AH1_2nd_14:bin_000063 p__Bacillota_A c__Clostridia o__Lachnospirales f__Lachnospiraceae g__Blautia_A s__
5 AH1_2nd_15:bin_000001 p__Pseudomonadota c__Alphaproteobacteria o__Rhizobiales f__Rhizobiaceae g__Agrobacterium s__Agrobacte…
6 AH1_2nd_16:bin_000033 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__ g__ s__
7 AH1_2nd_16:bin_000096 p__Bacillota_A c__Clostridia o__Oscillospirales f__Ruminococcaceae g__Ruthenibacterium s__
8 AH1_2nd_17:bin_000015 p__Bacillota c__Bacilli o__Erysipelotrichales f__Erysipelotrichaceae g__ s__
9 AH1_2nd_17:bin_000023 p__Bacillota c__Bacilli o__Erysipelotrichales f__Coprobacillaceae g__Coprobacillus s__
10 AH1_2nd_18:bin_000039 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Rikenellaceae g__Alistipes s__
# ℹ 22 more rowsCommunity elements differences in structural zeros
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>4.2.5 Differential abundance analysis: composition
phylo_samples_warm <- sample_metadata %>%
filter(Population == "Hot_dry" & time_point %in% c("Wild", "Acclimation") )%>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome_warm <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples_warm)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy_warm <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome_warm)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered_warm <- phyloseq(phylo_genome_warm, phylo_taxonomy_warm, phylo_samples_warm)ancom_rand_output_warm = ancombc2(data = physeq_genome_filtered_warm,
assay_name = "counts",
tax_level = NULL,
fix_formula = "time_point",
rand_formula = "(1|individual)",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered_warm@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output_warm$res %>%
dplyr::select(taxon, lfc_time_pointWild, p_time_pointWild) %>%
filter(p_time_pointWild < 0.05) %>%
dplyr::arrange(p_time_pointWild) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_time_pointWild)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_pointWild, y=forcats::fct_reorder(genome,lfc_time_pointWild), fill=phylum)) + #forcats::fct_rev()
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
# coord_flip()+
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))
Phyla of the significant MAGs
ancombc_rand_table_mag%>%
filter(lfc_time_pointWild<0) %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) phylum sample_count
1 Bacteroidota 4
2 Bacillota 2
3 Bacillota_A 2
4 Cyanobacteriota 2
5 Desulfobacterota 1
6 Verrucomicrobiota 1 phylum genus species
1 Bacillota_A
2 Verrucomicrobiota Akkermansia
3 Bacteroidota Odoribacter
4 Bacteroidota Alistipes
5 Bacillota C-19
6 Cyanobacteriota Scatousia
7 Cyanobacteriota Scatousia
8 Bacteroidota Phocaeicola
9 Bacillota Coprobacillus
10 Desulfobacterota Lawsonia
11 Bacteroidota Egerieousia
12 Bacillota_A IOR16 ancombc_rand_table_mag%>%
filter(lfc_time_pointWild<0) %>%
count(genus, name = "sample_count") %>%
arrange(desc(sample_count)) genus sample_count
1 Scatousia 2
2 1
3 Akkermansia 1
4 Alistipes 1
5 C-19 1
6 Coprobacillus 1
7 Egerieousia 1
8 IOR16 1
9 Lawsonia 1
10 Odoribacter 1
11 Phocaeicola 14.2.6 Differences in the functional capacity
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
filter(Population == "Hot_dry" & treatment %in% c("Acclimation", "Wild"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.2.6.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.355 0.0228
2 Wild 0.317 0.0335MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94062, p-value = 0.297
Wilcoxon rank sum exact test
data: value by treatment
W = 68, p-value = 0.01419
alternative hypothesis: true location shift is not equal to 0element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-Wild)%>%
mutate(group_color = ifelse(Difference <0, "Wild","Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=c("#d57d2c50","#d57d2c")) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")
4.2.6.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.348 0.0158
2 Wild 0.327 0.0245MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.95485, p-value = 0.5061
Wilcoxon rank sum exact test
data: value by treatment
W = 61, p-value = 0.07701
alternative hypothesis: true location shift is not equal to 0function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))| Code_function | Acclimation | Wild | Function | Difference | treatment |
|---|---|---|---|---|---|
| S02 | 0.1699218 | 0.2472264 | Appendages | -0.0773046 | Wild |
4.3 Does the antibiotic administration change the microbial community?
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(time_point %in% c("Acclimation", "Antibiotics") ) 4.3.1 Alpha diversity
label_map <- c(
"Cold_wet" = "Cold population",
"Hot_dry" = "Warm-population",
"richness" = "Species Richness",
"neutral" = "Neutral Diversity",
"phylogenetic" = "Phylogenetic Diversity")
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
filter(time_point %in% c("Acclimation", "Antibiotics") ) %>%
ggplot(aes(y = value, x = time_point, color=Population, fill=Population)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
facet_grid(metric ~ Population, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ time_point*Population+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(5.1645) ( log )
Formula: richness ~ time_point * Population + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
434.8 446.1 -211.4 422.8 43
Scaled residuals:
Min 1Q Median 3Q Max
-1.61979 -0.56608 0.06074 0.51882 2.39923
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0.2193 0.4683
Number of obs: 49, groups: individual, 27
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.7459 0.1623 23.073 < 2e-16 ***
time_pointAntibiotics -1.1523 0.1862 -6.190 6.02e-10 ***
PopulationHot_dry 0.7055 0.2719 2.595 0.00946 **
time_pointAntibiotics:PopulationHot_dry -0.3586 0.3037 -1.181 0.23769
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) tm_pnA PpltH_
tm_pntAntbt -0.466
PpltnHt_dry -0.598 0.278
tm_pntA:PH_ 0.292 -0.611 -0.461$emmeans
time_point emmean SE df asymp.LCL asymp.UCL
Acclimation 4.10 0.136 Inf 3.83 4.36
Antibiotics 2.77 0.151 Inf 2.47 3.06
Results are averaged over the levels of: Population
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Acclimation - Antibiotics 1.33 0.152 Inf 8.751 <.0001
Results are averaged over the levels of: Population
Results are given on the log (not the response) scale. $emmeans
Population emmean SE df asymp.LCL asymp.UCL
Cold_wet 3.17 0.145 Inf 2.89 3.45
Hot_dry 3.70 0.196 Inf 3.31 4.08
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
Cold_wet - Hot_dry -0.526 0.243 Inf -2.168 0.0302
Results are averaged over the levels of: time_point
Results are given on the log (not the response) scale. Neutral
Model_neutral <- nlme::lme(fixed = neutral ~ time_point*Population, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
342.5018 353.3418 -165.2509
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 3.282619 7.923552
Fixed effects: neutral ~ time_point * Population
Value Std.Error DF t-value p-value
(Intercept) 19.11713 2.078908 25 9.195758 0e+00
time_pointAntibiotics -11.46317 2.832826 20 -4.046551 6e-04
PopulationHot_dry 24.97345 3.534826 25 7.064972 0e+00
time_pointAntibiotics:PopulationHot_dry -20.91701 4.793363 20 -4.363745 3e-04
Correlation:
(Intr) tm_pnA PpltH_
time_pointAntibiotics -0.633
PopulationHot_dry -0.588 0.373
time_pointAntibiotics:PopulationHot_dry 0.374 -0.591 -0.632
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.8258256 -0.5685215 -0.1867832 0.5823712 2.0651455
Number of Observations: 49
Number of Groups: 27 $emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 31.60 1.77 25 27.96 35.2
Antibiotics 9.68 1.87 20 5.77 13.6
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Antibiotics 21.9 2.4 20 9.147 <.0001
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment $emmeans
Population emmean SE df lower.CL upper.CL
Cold_wet 13.4 1.61 20 10.0 16.7
Hot_dry 27.9 2.22 20 23.3 32.5
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Cold_wet - Hot_dry -14.5 2.74 20 -5.288 <.0001
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment Phylogenetic
Model_phylo <- nlme::lme(fixed = phylogenetic ~ time_point*Population, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
203.6546 214.4946 -95.8273
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.7173285 1.688356
Fixed effects: phylogenetic ~ time_point * Population
Value Std.Error DF t-value p-value
(Intercept) 5.365728 0.4446271 25 12.067930 0.0000
time_pointAntibiotics -0.761455 0.6038652 20 -1.260969 0.2218
PopulationHot_dry 1.097291 0.7560383 25 1.451370 0.1591
time_pointAntibiotics:PopulationHot_dry -0.702255 1.0216421 20 -0.687379 0.4997
Correlation:
(Intr) tm_pnA PpltH_
time_pointAntibiotics -0.631
PopulationHot_dry -0.588 0.371
time_pointAntibiotics:PopulationHot_dry 0.373 -0.591 -0.629
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.85993684 -0.51436620 -0.05130501 0.48944142 1.80209573
Number of Observations: 49
Number of Groups: 27 $emmeans
time_point emmean SE df lower.CL upper.CL
Acclimation 31.60 1.77 25 27.96 35.2
Antibiotics 9.68 1.87 20 5.77 13.6
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Antibiotics 21.9 2.4 20 9.147 <.0001
Results are averaged over the levels of: Population
Degrees-of-freedom method: containment $emmeans
Population emmean SE df lower.CL upper.CL
Cold_wet 13.4 1.61 20 10.0 16.7
Hot_dry 27.9 2.22 20 23.3 32.5
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Cold_wet - Hot_dry -14.5 2.74 20 -5.288 <.0001
Results are averaged over the levels of: time_point
Degrees-of-freedom method: containment 4.3.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(time_point %in% c("Acclimation", "Antibiotics")) %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(time_point %in% c("Acclimation", "Antibiotics"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.03530 0.035300 10.073 999 0.004 **
Residuals 47 0.16471 0.003505
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(richness ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 1.9348458 | 0.1002347 | 6.065418 | 0.001 |
| Population | 1 | 2.1356198 | 0.1106358 | 6.694811 | 0.001 |
| time_point:Population | 1 | 0.8778553 | 0.0454773 | 2.751930 | 0.002 |
| Residual | 45 | 14.3548318 | 0.7436522 | NA | NA |
| Total | 48 | 19.3031527 | 1.0000000 | NA | NA |
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(richness),]
pairwise <- pairwise.adonis(richness, subset_meta_arrange$Population, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_wet vs Hot_dry | 1 | 2.135115 | 5.845187 | 0.1106096 | 0.001 | 0.001 | ** |
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Acclimation vs Antibiotics | 1 | 1.934846 | 5.235844 | 0.1002347 | 0.001 | 0.001 | ** |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.051657 0.051657 9.9224 999 0.003 **
Residuals 47 0.244689 0.005206
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(neutral ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 2.0429592 | 0.11061666 | 7.254623 | 0.001 |
| Population | 1 | 2.8764490 | 0.15574623 | 10.214376 | 0.001 |
| time_point:Population | 1 | 0.8770558 | 0.04748846 | 3.114457 | 0.001 |
| Residual | 45 | 12.6723552 | 0.68614864 | NA | NA |
| Total | 48 | 18.4688192 | 1.00000000 | NA | NA |
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(neutral),]
pairwise <- pairwise.adonis(neutral, subset_meta_arrange$Population, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_wet vs Hot_dry | 1 | 2.876596 | 8.670991 | 0.1557542 | 0.001 | 0.001 | ** |
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Acclimation vs Antibiotics | 1 | 2.042959 | 5.845605 | 0.1106167 | 0.001 | 0.001 | ** |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.67439 0.67439 55.932 999 0.001 ***
Residuals 47 0.56670 0.01206
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(phylo ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 1.8783556 | 0.23032747 | 16.358708 | 0.001 |
| Population | 1 | 0.7551942 | 0.09260332 | 6.577030 | 0.001 |
| time_point:Population | 1 | 0.3545686 | 0.04347786 | 3.087959 | 0.022 |
| Residual | 45 | 5.1670341 | 0.63359135 | NA | NA |
| Total | 48 | 8.1551525 | 1.00000000 | NA | NA |
subset_meta_arrange <- column_to_rownames(subset_meta, "Tube_code")
subset_meta_arrange<-subset_meta_arrange[labels(phylo),]
pairwise <- pairwise.adonis(phylo, subset_meta_arrange$Population, perm=999)
pairwise%>%
tt()| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Cold_wet vs Hot_dry | 1 | 0.753589 | 4.785297 | 0.09240649 | 0.001 | 0.001 | ** |
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| Acclimation vs Antibiotics | 1 | 1.878356 | 14.06493 | 0.2303275 | 0.001 | 0.001 | ** |
dbRDA
#Richness
cca_ord <- capscale(formula = richness ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta$time_pointAntibiotics" = "Antibiotics",
"subset_meta$PopulationHot_dry"="Warm population",
"subset_meta$time_pointAntibiotics:subset_meta$PopulationHot_dry"="Interaction Wam population")
beta_richness_rda <-CAP_df %>%
group_by(Population, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Neutral
cca_ord <- capscale(formula = neutral ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta$time_pointAntibiotics" = "Antibiotics",
"subset_meta$PopulationHot_dry"="Warm population",
"subset_meta$time_pointAntibiotics:subset_meta$PopulationHot_dry"="Interaction Wam population")
beta_neutral_rda <- CAP_df %>%
group_by(Population, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50"))+
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()
#Phylogenetic
cca_ord <- capscale(formula = phylo ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(biplot_scores$Variable,
"subset_meta$time_pointAntibiotics" = "Antibiotics",
"subset_meta$PopulationHot_dry"="Warm population",
"subset_meta$time_pointAntibiotics:subset_meta$PopulationHot_dry"="Interaction Wam population")
beta_phylogenetic_rda<- CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()ggarrange(beta_richness_rda, beta_neutral_rda, beta_phylogenetic_rda, ncol=3, nrow=1, common.legend = TRUE, legend="right")
4.3.3 Cold population
4.3.3.1 Structural zeros
acclimation_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point == "Acclimation") %>%
dplyr::select(sample) %>%
pull()
antibiotics_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point == "Antibiotics") %>%
dplyr::select(sample) %>% pull()
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",subset_meta$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>%
mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>%
mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>%
mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>%
filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE) %>%
mutate(present = case_when(
all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
!all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
!all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Acclimation
structural_zeros %>%
filter(present=="acclimation") %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) # A tibble: 9 × 2
# Rowwise:
phylum sample_count
<chr> <int>
1 p__Bacillota_A 43
2 p__Bacteroidota 15
3 p__Bacillota 8
4 p__Pseudomonadota 7
5 p__Cyanobacteriota 3
6 p__Verrucomicrobiota 2
7 p__Bacillota_B 1
8 p__Bacillota_C 1
9 p__Fusobacteriota 1Antibiotics
# A tibble: 4 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_7:bin_000003 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Proteus s__Proteus cibarius
2 LI1_2nd_7:bin_000001 p__Bacillota_A c__Clostridia o__Clostridiales f__Clostridiaceae g__Sarcina s__
3 AH1_2nd_7:bin_000055 p__Bacillota c__Bacilli o__Mycoplasmatales f__Mycoplasmoidaceae g__Ureaplasma s__
4 AH1_2nd_13:bin_000025 p__Bacillota_A c__Clostridia o__Christensenellales f__UBA3700 g__ s__ Community elements differences in structural zeros
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>4.3.3.2 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point %in% c("Acclimation", "Antibiotics") )%>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "time_point",
rand_formula = "(1|individual)",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output$res %>%
dplyr::select(taxon, lfc_time_pointAntibiotics, p_time_pointAntibiotics) %>%
filter(p_time_pointAntibiotics < 0.05) %>%
dplyr::arrange(p_time_pointAntibiotics) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_time_pointAntibiotics)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_pointAntibiotics, y=forcats::fct_reorder(genome,lfc_time_pointAntibiotics), fill=phylum)) + #forcats::fct_rev()
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
# coord_flip()+
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))
Phyla of the significant MAGs
ancombc_rand_table_mag%>%
filter(lfc_time_pointAntibiotics<0) %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) phylum sample_count
1 Bacteroidota 15
2 Bacillota_A 5
3 Bacillota 1
4 Campylobacterota 1
5 Cyanobacteriota 1 phylum genus species
1 Campylobacterota Helicobacter_J
2 Bacillota_A
3 Bacteroidota Odoribacter
4 Bacteroidota Bacteroides
5 Bacteroidota Bacteroides
6 Bacillota
7 Bacteroidota Bacteroides
8 Bacteroidota Phocaeicola
9 Bacteroidota Odoribacter
10 Bacteroidota Phocaeicola
11 Bacteroidota Bacteroides
12 Bacteroidota Bacteroides
13 Bacteroidota Parabacteroides
14 Bacteroidota
15 Bacteroidota Parabacteroides
16 Cyanobacteriota Scatousia
17 Bacillota_A
18 Bacillota_A Lacrimispora
19 Bacteroidota Alistipes
20 Bacillota_A
21 Bacillota_A
22 Bacteroidota Alistipes
23 Bacteroidota Odoribacter ancombc_rand_table_mag%>%
filter(lfc_time_pointAntibiotics<0) %>%
count(genus, name = "sample_count") %>%
arrange(desc(sample_count)) genus sample_count
1 6
2 Bacteroides 5
3 Odoribacter 3
4 Alistipes 2
5 Parabacteroides 2
6 Phocaeicola 2
7 Helicobacter_J 1
8 Lacrimispora 1
9 Scatousia 14.3.3.3 Differences in the functional capacity
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Antibiotics"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.3.3.3.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.335 0.0324
2 Antibiotics 0.247 0.136 MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94332, p-value = 0.09304
Wilcoxon rank sum exact test
data: value by treatment
W = 190, p-value = 0.01768
alternative hypothesis: true location shift is not equal to 0element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-Antibiotics)%>%
mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=c('#008080', '#00808050')) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")
4.3.3.3.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.330 0.0320
2 Antibiotics 0.256 0.126 MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.95742, p-value = 0.2331
Wilcoxon rank sum exact test
data: value by treatment
W = 188, p-value = 0.02195
alternative hypothesis: true location shift is not equal to 0function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))| Code_function | Acclimation | Antibiotics | Function | Difference | treatment |
|---|---|---|---|---|---|
| B06 | 0.68110280 | 0.47325490 | Organic anion biosynthesis | 0.20784790 | Acclimation |
| B02 | 0.59963040 | 0.41562100 | Amino acid biosynthesis | 0.18400940 | Acclimation |
| D02 | 0.38587770 | 0.20636760 | Polysaccharide degradation | 0.17951010 | Acclimation |
| S03 | 0.27103471 | 0.09357224 | Spore | 0.17746247 | Acclimation |
| B01 | 0.84215140 | 0.69078910 | Nucleic acid biosynthesis | 0.15136230 | Acclimation |
| B07 | 0.44558700 | 0.30432910 | Vitamin biosynthesis | 0.14125790 | Acclimation |
| B08 | 0.44596150 | 0.32044710 | Aromatic compound biosynthesis | 0.12551440 | Acclimation |
| D09 | 0.25533360 | 0.13438350 | Antibiotic degradation | 0.12095010 | Acclimation |
| B04 | 0.54457570 | 0.42921780 | SCFA biosynthesis | 0.11535790 | Acclimation |
| D03 | 0.29210750 | 0.20698550 | Sugar degradation | 0.08512200 | Acclimation |
| D05 | 0.28856110 | 0.22258820 | Amino acid degradation | 0.06597290 | Acclimation |
| B10 | 0.05947806 | 0.03621687 | Antibiotic biosynthesis | 0.02326119 | Acclimation |
4.3.4 Warm population
4.3.4.1 Structural zeros
acclimation_samples <- sample_metadata %>%
filter(Population == "Hot_dry" & time_point == "Acclimation") %>%
dplyr::select(sample) %>%
pull()
antibiotics_samples <- sample_metadata %>%
filter(Population == "Hot_dry" & time_point == "Antibiotics") %>%
dplyr::select(sample) %>% pull()
sample_sub <- sample_metadata %>%
filter(Population == "Hot_dry" & treatment %in% c("Acclimation", "Antibiotics"))
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",sample_sub$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>%
mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>%
mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>%
mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>%
filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE) %>%
mutate(present = case_when(
all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
!all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
!all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Acclimation
structural_zeros %>%
filter(present=="acclimation") %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) # A tibble: 12 × 2
# Rowwise:
phylum sample_count
<chr> <int>
1 p__Bacillota_A 63
2 p__Bacteroidota 32
3 p__Bacillota 16
4 p__Cyanobacteriota 9
5 p__Pseudomonadota 7
6 p__Bacillota_B 3
7 p__Desulfobacterota 3
8 p__Bacillota_C 2
9 p__Verrucomicrobiota 2
10 p__Actinomycetota 1
11 p__Campylobacterota 1
12 p__Elusimicrobiota 1Antibiotics
# A tibble: 7 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_7:bin_000003 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Proteus s__Proteus cibarius
2 LI1_2nd_8:bin_000030 p__Actinomycetota c__Actinomycetia o__Mycobacteriales f__Mycobacteriaceae g__Corynebacterium s__
3 LI1_2nd_8:bin_000077 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Tannerellaceae g__Parabacteroides s__Parabacteroides…
4 LI1_2nd_5:bin_000032 p__Bacillota_A c__Clostridia o__Christensenellales f__UBA1242 g__Caccovivens s__
5 AH1_2nd_18:bin_000013 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Tannerellaceae g__Parabacteroides s__
6 LI1_2nd_5:bin_000069 p__Bacillota c__Bacilli o__RF39 f__UBA660 g__ s__
7 AH1_2nd_20:bin_000061 p__Bacillota c__Bacilli o__Lactobacillales f__Enterococcaceae g__Enterococcus s__ Community elements differences in structural zeros
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)# A tibble: 20 × 3
trait p_value p_adjust
<chr> <dbl> <dbl>
1 B0220 0.00618 0.0497
2 B0708 0.00280 0.0497
3 B0709 0.00662 0.0497
4 B1012 0.00618 0.0497
5 D0203 0.00559 0.0497
6 D0205 0.00280 0.0497
7 D0206 0.00280 0.0497
8 D0207 0.00280 0.0497
9 D0210 0.00280 0.0497
10 D0211 0.00662 0.0497
11 D0213 0.00685 0.0497
12 D0305 0.00685 0.0497
13 D0308 0.00685 0.0497
14 D0610 0.00618 0.0497
15 D0706 0.00662 0.0497
16 D0902 0.00618 0.0497
17 D0905 0.00618 0.0497
18 D0908 0.00280 0.0497
19 D0911 0.00618 0.0497
20 S0104 0.00685 0.04974.3.4.2 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
filter(Population == "Hot_dry" & time_point %in% c("Acclimation", "Antibiotics") )%>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "time_point",
rand_formula = "(1|individual)",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output$res %>%
dplyr::select(taxon, lfc_time_pointAntibiotics, p_time_pointAntibiotics) %>%
filter(p_time_pointAntibiotics < 0.05) %>%
dplyr::arrange(p_time_pointAntibiotics) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_time_pointAntibiotics)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_pointAntibiotics, y=forcats::fct_reorder(genome,lfc_time_pointAntibiotics), fill=phylum)) +
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))Phyla of the significant MAGs
ancombc_rand_table_mag%>%
filter(lfc_time_pointAntibiotics<0) %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count))
ancombc_rand_table_mag%>%
filter(lfc_time_pointAntibiotics<0) %>%
select(phylum, genus, species)
ancombc_rand_table_mag%>%
filter(lfc_time_pointAntibiotics<0) %>%
count(genus, name = "sample_count") %>%
arrange(desc(sample_count))4.3.4.3 Differences in the functional capacity
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Antibiotics"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.3.4.3.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.335 0.0324
2 Antibiotics 0.247 0.136 MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94332, p-value = 0.09304
Wilcoxon rank sum exact test
data: value by treatment
W = 190, p-value = 0.01768
alternative hypothesis: true location shift is not equal to 0element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-Antibiotics)%>%
mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=c("#d57d2c50","#d57d2c")) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")
4.3.4.3.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.330 0.0320
2 Antibiotics 0.256 0.126 MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.95742, p-value = 0.2331
Wilcoxon rank sum exact test
data: value by treatment
W = 188, p-value = 0.02195
alternative hypothesis: true location shift is not equal to 0function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))| Code_function | Acclimation | Antibiotics | Function | Difference | treatment |
|---|---|---|---|---|---|
| B06 | 0.68110280 | 0.47325490 | Organic anion biosynthesis | 0.20784790 | Acclimation |
| B02 | 0.59963040 | 0.41562100 | Amino acid biosynthesis | 0.18400940 | Acclimation |
| D02 | 0.38587770 | 0.20636760 | Polysaccharide degradation | 0.17951010 | Acclimation |
| S03 | 0.27103471 | 0.09357224 | Spore | 0.17746247 | Acclimation |
| B01 | 0.84215140 | 0.69078910 | Nucleic acid biosynthesis | 0.15136230 | Acclimation |
| B07 | 0.44558700 | 0.30432910 | Vitamin biosynthesis | 0.14125790 | Acclimation |
| B08 | 0.44596150 | 0.32044710 | Aromatic compound biosynthesis | 0.12551440 | Acclimation |
| D09 | 0.25533360 | 0.13438350 | Antibiotic degradation | 0.12095010 | Acclimation |
| B04 | 0.54457570 | 0.42921780 | SCFA biosynthesis | 0.11535790 | Acclimation |
| D03 | 0.29210750 | 0.20698550 | Sugar degradation | 0.08512200 | Acclimation |
| D05 | 0.28856110 | 0.22258820 | Amino acid degradation | 0.06597290 | Acclimation |
| B10 | 0.05947806 | 0.03621687 | Antibiotic biosynthesis | 0.02326119 | Acclimation |
4.4 Does antibiotic administration remove the differences found in the two populations?
4.4.1 Shapiro test
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "Antibiotics" ) %>%
filter(metric=="richness") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.001165318alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "Antibiotics" ) %>%
filter(metric=="neutral") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.0003462674alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "Antibiotics" ) %>%
filter(metric=="phylogenetic") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.066284424.4.2 Alpha diversity
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(time_point == "Antibiotics" )alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness", "neutral"))) %>%
filter(metric!="phylogenetic") %>%
filter(time_point == "Antibiotics" ) %>%
ggplot(aes(y = value, x = Population, color=Population, fill=Population)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
stat_compare_means(method = "t.test", show.legend = F, size = 3, label.x = c(1.5))+
facet_grid( ~ metric, scales = "free_y")+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral"))) %>%
filter(metric!="phylogenetic") %>%
filter(time_point == "Antibiotics" ) %>%
ggplot(aes(y = value, x = Population, color=Population, fill=Population)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
4.4.3 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(time_point == "Antibiotics") %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(time_point == "Antibiotics")Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$Population) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.015377 0.0153774 6.8934 999 0.02 *
Residuals 21 0.046845 0.0022307
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(richness ~ Population,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.361743 | 0.1533845 | 3.80465 | 0.001 |
| Residual | 21 | 7.516224 | 0.8466155 | NA | NA |
| Total | 22 | 8.877966 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$Population) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.030649 0.0306488 3.8694 999 0.061 .
Residuals 21 0.166339 0.0079209
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(neutral ~ Population,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.787698 | 0.208772 | 5.541022 | 0.001 |
| Residual | 21 | 6.775221 | 0.791228 | NA | NA |
| Total | 22 | 8.562919 | 1.000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$Population) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.012165 0.012165 0.9979 999 0.328
Residuals 21 0.256012 0.012191 adonis2(phylo ~ Population,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.8970751 | 0.1890846 | 4.896659 | 0.001 |
| Residual | 21 | 3.8472307 | 0.8109154 | NA | NA |
| Total | 22 | 4.7443057 | 1.0000000 | NA | NA |
4.5 Are the microbial communities similar in both donor samples?
4.5.1 Alpha diversity
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type =="Hot_control" & time_point %in% c( "Transplant1", "Transplant2"))samples_to_keep <- sample_metadata %>%
filter(type =="Hot_control" & time_point %in% c( "Transplant1", "Transplant2"))%>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type =="Hot_control" & time_point %in% c( "Transplant1", "Transplant2"))alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
filter(type =="Hot_control" & time_point %in% c( "Transplant1", "Transplant2")) %>%
ggplot(aes(y = value, x = time_point, color=time_point, fill=time_point)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="time_point",
breaks=c("Transplant1","Transplant2"),
values=c('#d5992c', "#d5b52c")) +
scale_fill_manual(name="time_point",
breaks=c("Transplant1","Transplant2"),
values=c('#d5992c50', "#d5b52c50")) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(17.9668) ( log )
Formula: richness ~ time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
150.4 153.3 -71.2 142.4 11
Scaled residuals:
Min 1Q Median 3Q Max
-2.02956 -0.52933 0.03513 0.56378 1.56130
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0.009989 0.09995
Number of obs: 15, groups: individual, 8
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.60703 0.09905 46.514 <2e-16 ***
time_pointTransplant2 0.05482 0.13299 0.412 0.68
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
tm_pntTrns2 -0.624Neutral
Model_neutral <- nlme::lme(fixed = neutral ~ time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
115.9183 118.1781 -53.95914
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 9.533528 10.15744
Fixed effects: neutral ~ time_point
Value Std.Error DF t-value p-value
(Intercept) 43.94053 4.925212 7 8.921551 0.0000
time_pointTransplant2 4.31623 5.338413 6 0.808523 0.4497
Correlation:
(Intr)
time_pointTransplant2 -0.491
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.57714267 -0.56764018 0.05406399 0.55332481 1.12750551
Number of Observations: 15
Number of Groups: 8 Phylogenetic
Model_phylo <- nlme::lme(fixed = phylogenetic ~ time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
41.16537 43.42517 -16.58269
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 1.16159 0.3042944
Fixed effects: phylogenetic ~ time_point
Value Std.Error DF t-value p-value
(Intercept) 5.785003 0.4245418 7 13.626462 0.0000
time_pointTransplant2 -0.018098 0.1623255 6 -0.111494 0.9149
Correlation:
(Intr)
time_pointTransplant2 -0.168
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.2428123 -0.3756145 -0.0502414 0.4569146 1.3175579
Number of Observations: 15
Number of Groups: 8 4.5.2 Beta diversity
Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00014 0.0001396 0.0173 999 0.906
Residuals 13 0.10481 0.0080621 adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.06889268 | 0.03010712 | 0.4035421 | 0.9375 |
| Residual | 13 | 2.21935895 | 0.96989288 | NA | NA |
| Total | 14 | 2.28825162 | 1.00000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000655 0.0006546 0.0514 999 0.836
Residuals 13 0.165409 0.0127237 adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.09294671 | 0.03882177 | 0.5250671 | 0.7265625 |
| Residual | 13 | 2.30124351 | 0.96117823 | NA | NA |
| Total | 14 | 2.39419022 | 1.00000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.003691 0.0036912 0.3071 999 0.673
Residuals 13 0.156266 0.0120205 adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.009773632 | 0.01921019 | 0.2546239 | 0.7734375 |
| Residual | 13 | 0.498999565 | 0.98078981 | NA | NA |
| Total | 14 | 0.508773196 | 1.00000000 | NA | NA |
4.6 Does the donor sample maintain the microbial community found in acclimation?
sample_metadata <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Transplant1", "Transplant2") ~ "Transplant",
TRUE ~ treatment
))
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant"))sample_metadata <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Transplant1", "Transplant2") ~ "Transplant",
TRUE ~ treatment
))
samples_to_keep <- sample_metadata %>%
filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant"))4.6.1 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c("#d57d2c", "#d5b52c")) +
scale_fill_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c("#d57d2c50", "#d5b52c50")) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(19.6072) ( log )
Formula: richness ~ treatment + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
228.0 232.7 -110.0 220.0 20
Scaled residuals:
Min 1Q Median 3Q Max
-2.37872 -0.53180 0.06467 0.46904 1.76837
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0 0
Number of obs: 24, groups: individual, 9
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.4569 0.0834 53.441 <2e-16 ***
treatmentTransplant 0.1855 0.1049 1.769 0.0769 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
trtmntTrnsp -0.795
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Neutral
Model_neutral <- nlme::lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
185.3367 189.7009 -88.66837
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.001935042 12.18199
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 44.09058 4.060663 14 10.857976 0.0000
treatmentTransplant 1.80350 5.136377 14 0.351122 0.7307
Correlation:
(Intr)
treatmentTransplant -0.791
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.66610957 -0.48314823 -0.03902931 0.62479437 2.02597977
Number of Observations: 24
Number of Groups: 9 Phylogenetic
Model_phylo <- nlme::lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
80.6646 85.02877 -36.3323
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.7197786 0.9593315
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 6.463020 0.3997775 14 16.166543 0.0000
treatmentTransplant -0.577492 0.4112749 14 -1.404151 0.1821
Correlation:
(Intr)
treatmentTransplant -0.622
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.329685757 -0.568101033 -0.001989183 0.552370188 1.529548996
Number of Observations: 24
Number of Groups: 9 4.6.2 Beta diversity
Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00836 0.0083600 1.4409 999 0.274
Residuals 22 0.12764 0.0058019 adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.2657351 | 0.06396199 | 1.503319 | 0.102 |
| Residual | 22 | 3.8888430 | 0.93603801 | NA | NA |
| Total | 23 | 4.1545781 | 1.00000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000661 0.0006614 0.0736 999 0.798
Residuals 22 0.197804 0.0089911 adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3164816 | 0.07596593 | 1.808646 | 0.083 |
| Residual | 22 | 3.8496178 | 0.92403407 | NA | NA |
| Total | 23 | 4.1660995 | 1.00000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00204 0.0020405 0.2622 999 0.685
Residuals 22 0.17118 0.0077807 adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.0412056 | 0.056244 | 1.31111 | 0.25 |
| Residual | 22 | 0.6914166 | 0.943756 | NA | NA |
| Total | 23 | 0.7326222 | 1.000000 | NA | NA |
p0<-richness %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c("#d57d2c", "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p1<-neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c("#d57d2c", "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p2<-phylo %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c("#d57d2c", "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
4.7 Is the donor sample microbiota different to recipient microbiota?
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant"))4.7.1 Alpha diversity
sample_metadata <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Transplant1", "Transplant2") ~ "Transplant",
TRUE ~ treatment
))
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Treatment" & treatment %in% c("Acclimation","Transplant"))alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c('#008080', "#d5b52c")) +
scale_fill_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c('#00808050', "#d5b52c50")) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(6.4835) ( log )
Formula: richness ~ treatment + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
219.3 223.6 -105.6 211.3 18
Scaled residuals:
Min 1Q Median 3Q Max
-1.9573 -0.5447 0.1304 0.5773 1.2167
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 5.302e-13 7.281e-07
Number of obs: 22, groups: individual, 9
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.8067 0.1400 27.19 < 2e-16 ***
treatmentTransplant 0.8399 0.1795 4.68 2.87e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
trtmntTrnsp -0.780
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Neutral
Model_neutral <-nlme:: lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
169.583 173.5659 -80.79148
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 4.713488 11.40939
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 17.31015 4.114893 12 4.206708 0.0012
treatmentTransplant 28.70744 5.016259 12 5.722879 0.0001
Correlation:
(Intr)
treatmentTransplant -0.701
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.4861810 -0.5958348 0.1169620 0.5490436 1.7889456
Number of Observations: 22
Number of Groups: 9 Phylogenetic
Model_phylo <- nlme::lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
80.93175 84.91468 -36.46587
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.8689563 1.118803
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 5.488256 0.4722057 12 11.622595 0.00
treatmentTransplant 0.226402 0.5020134 12 0.450988 0.66
Correlation:
(Intr)
treatmentTransplant -0.587
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.42844864 -0.44112019 -0.05260838 0.51238319 1.55109108
Number of Observations: 22
Number of Groups: 9 4.7.2 Beta diversity
Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.12795 0.127946 13.179 999 0.001 ***
Residuals 20 0.19416 0.009708
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.777815 | 0.2760196 | 7.625059 | 0.002 |
| Residual | 20 | 4.663086 | 0.7239804 | NA | NA |
| Total | 21 | 6.440902 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.071502 0.071502 5.4967 999 0.032 *
Residuals 20 0.260166 0.013008
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 2.112572 | 0.3210813 | 9.458609 | 0.005 |
| Residual | 20 | 4.466983 | 0.6789187 | NA | NA |
| Total | 21 | 6.579556 | 1.0000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04731 0.047307 2.6157 999 0.116
Residuals 20 0.36172 0.018086 adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3778248 | 0.239656 | 6.303884 | 0.002 |
| Residual | 20 | 1.1987049 | 0.760344 | NA | NA |
| Total | 21 | 1.5765298 | 1.000000 | NA | NA |
p0<-richness %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c('#008080', "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p1<-neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c('#008080', "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p2<-phylo %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation","Transplant"),
values=c('#008080', "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
4.8 Does FMT change the microbial community over time?
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Treatment" & treatment %in% c("FMT1","FMT2") ) 4.8.1 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral", "phylogenetic"))) %>%
filter(type=="Treatment" & treatment %in% c("FMT1", "FMT2")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="treatment",
values=c("#76b183", "#2b4b32")) +
scale_fill_manual(name="treatment",
values=c("#76b18350", "#2b4b3250")) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(4.3876) ( log )
Formula: richness ~ treatment + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
171.0 174.3 -81.5 163.0 13
Scaled residuals:
Min 1Q Median 3Q Max
-1.73495 -0.71265 -0.07086 0.40744 1.88240
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 1.073e-11 3.275e-06
Number of obs: 17, groups: individual, 9
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9343 0.1759 22.369 <2e-16 ***
treatmentFMT2 0.4052 0.2402 1.687 0.0917 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
tretmntFMT2 -0.732
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Neutral
Model_neutral <- nlme::lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
128.7946 131.6268 -60.3973
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 3.497472 11.25384
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 24.65947 4.164756 8 5.920988 0.0004
treatmentFMT2 14.87494 5.482531 7 2.713151 0.0301
Correlation:
(Intr)
treatmentFMT2 -0.7
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.76462324 -0.55701822 -0.04763166 0.55267588 1.30333436
Number of Observations: 17
Number of Groups: 9 Phylogenetic
Model_phylo <- nlme::lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
56.05281 58.88501 -24.0264
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.6954001 0.8350046
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 4.154090 0.3805932 8 10.914778 0.0000
treatmentFMT2 0.912624 0.4105974 7 2.222673 0.0616
Correlation:
(Intr)
treatmentFMT2 -0.583
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.1361565 -0.5779234 -0.1465750 0.4815529 2.3352020
Number of Observations: 17
Number of Groups: 9 4.8.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("FMT1", "FMT2")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("FMT1", "FMT2"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.017610 0.017610 2.8225 999 0.115
Residuals 15 0.093585 0.006239 adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3560084 | 0.07968734 | 1.298809 | 0.00390625 |
| Residual | 15 | 4.1115570 | 0.92031266 | NA | NA |
| Total | 16 | 4.4675655 | 1.00000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.009828 0.0098278 0.7113 999 0.441
Residuals 15 0.207263 0.0138175 adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3149954 | 0.08110541 | 1.323962 | 0.00390625 |
| Residual | 15 | 3.5687823 | 0.91889459 | NA | NA |
| Total | 16 | 3.8837776 | 1.00000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.010955 0.010955 0.6602 999 0.515
Residuals 15 0.248892 0.016593 adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.06391536 | 0.09449338 | 1.565312 | 0.0703125 |
| Residual | 15 | 0.61248501 | 0.90550662 | NA | NA |
| Total | 16 | 0.67640037 | 1.00000000 | NA | NA |
p0<-richness %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("FMT1", "FMT2"),
values=c("#76b183", "#2b4b32")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p1<-neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("FMT1", "FMT2"),
values=c("#76b183", "#2b4b32")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p2<-phylo %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("FMT1", "FMT2"),
values=c("#76b183", "#2b4b32")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
4.9 Do FMT change the microbial community compared to antibiotics and acclimation?
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Treatment" & treatment %in% c("Antibiotics", "Acclimation","FMT1") ) 4.9.1 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral", "phylogenetic"))) %>%
filter(type=="Treatment" & treatment %in% c( "Antibiotics","Acclimation", "FMT1")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(name="Treatment",
breaks=c("Antibiotics","Acclimation", "FMT1"),
values=c("#003a45",'#008080', "#76b183")) +
scale_fill_manual(name="Treatment",
breaks=c("Antibiotics","Acclimation", "FMT1"),
values=c("#003a4550",'#00808050',"#76b18350")) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness: Problems to calculate
Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Neutral
Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
173.0404 178.263 -81.5202
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 4.263247 9.322271
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 17.310151 3.416951 13 5.065963 0.0002
treatmentAntibiotics -8.868175 4.744329 13 -1.869216 0.0843
treatmentFMT1 7.289358 4.552796 13 1.601073 0.1334
Correlation:
(Intr) trtmnA
treatmentAntibiotics -0.596
treatmentFMT1 -0.621 0.457
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.93185401 -0.56384573 -0.04015247 0.47657219 1.55593253
Number of Observations: 24
Number of Groups: 9 Phylogenetic
Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
92.7946 98.01721 -41.3973
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.5507824 1.405346
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 5.488256 0.5031412 13 10.907984 0.0000
treatmentAntibiotics -1.198294 0.7137109 13 -1.678962 0.1170
treatmentFMT1 -1.351708 0.6855445 13 -1.971729 0.0703
Correlation:
(Intr) trtmnA
treatmentAntibiotics -0.611
treatmentFMT1 -0.636 0.456
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.00120146 -0.44426020 0.00785757 0.38719755 1.97397308
Number of Observations: 24
Number of Groups: 9 4.9.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Antibiotics", "FMT1")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Antibiotics", "FMT1"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.01582 0.0079101 1.05 999 0.398
Residuals 21 0.15821 0.0075336 adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.840841 | 0.2018508 | 2.655435 | 0.001 |
| Residual | 21 | 7.278967 | 0.7981492 | NA | NA |
| Total | 23 | 9.119808 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.8287871 2.293722 0.1407734 0.003 0.009 *
2 Acclimation vs FMT1 1 0.9572484 2.939042 0.1638350 0.001 0.003 *
3 Antibiotics vs FMT1 1 0.9764237 2.751191 0.1746656 0.005 0.015 .Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.01796 0.0089778 0.5959 999 0.535
Residuals 21 0.31640 0.0150667 adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 2.357142 | 0.2699241 | 3.882065 | 0.001 |
| Residual | 21 | 6.375470 | 0.7300759 | NA | NA |
| Total | 23 | 8.732613 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.8814932 2.747044 0.1640316 0.001 0.003 *
2 Acclimation vs FMT1 1 1.3133104 4.634354 0.2360329 0.001 0.003 *
3 Antibiotics vs FMT1 1 1.3427497 4.355528 0.2509591 0.001 0.003 *Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.15410 0.077052 2.6434 999 0.101
Residuals 21 0.61214 0.029149 adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.327076 | 0.3647564 | 6.029092 | 0.001 |
| Residual | 21 | 2.311178 | 0.6352436 | NA | NA |
| Total | 23 | 3.638254 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.6885926 5.124832 0.2679674 0.004 0.012 .
2 Acclimation vs FMT1 1 0.2548027 3.292748 0.1800029 0.038 0.114
3 Antibiotics vs FMT1 1 1.1000473 9.048069 0.4103792 0.004 0.012 .p0<-richness %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation", "Antibiotics", "FMT1"),
values=c("#003a45",'#008080', "#76b183")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p1<-neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation", "Antibiotics", "FMT1"),
values=c("#003a45",'#008080', "#76b183")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p2<-phylo %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("Acclimation", "Antibiotics", "FMT1"),
values=c("#003a45",'#008080', "#76b183")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
4.10 Is the gut microbiota similar to the donor after FMT?
sample_metadata <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("FMT1", "FMT2") ~ "FMT",
TRUE ~ treatment
))
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))4.10.1 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("neutral","richness","phylogenetic"))) %>%
filter(type=="Treatment" & treatment %in% c( "Transplant", "FMT")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(name="treatment",
breaks=c("FMT","Transplant"),
values=c('#40714b', "#d5b52c")) +
scale_fill_manual(name="treatment",
breaks=c("FMT","Transplant"),
values=c('#40714b50', "#d5b52c50")) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
Richness
Modelq0GLMMNB <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ treatment)Neutral
Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
173.0404 178.263 -81.5202
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 4.263247 9.322271
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 17.310151 3.416951 13 5.065963 0.0002
treatmentAntibiotics -8.868175 4.744329 13 -1.869216 0.0843
treatmentFMT1 7.289358 4.552796 13 1.601073 0.1334
Correlation:
(Intr) trtmnA
treatmentAntibiotics -0.596
treatmentFMT1 -0.621 0.457
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.93185401 -0.56384573 -0.04015247 0.47657219 1.55593253
Number of Observations: 24
Number of Groups: 9 Phylogenetic
Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
92.7946 98.01721 -41.3973
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.5507824 1.405346
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 5.488256 0.5031412 13 10.907984 0.0000
treatmentAntibiotics -1.198294 0.7137109 13 -1.678962 0.1170
treatmentFMT1 -1.351708 0.6855445 13 -1.971729 0.0703
Correlation:
(Intr) trtmnA
treatmentAntibiotics -0.611
treatmentFMT1 -0.636 0.456
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.00120146 -0.44426020 0.00785757 0.38719755 1.97397308
Number of Observations: 24
Number of Groups: 9 4.10.2 Beta diversity
Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.11100 0.111002 11.752 999 0.001 ***
Residuals 28 0.26447 0.009445
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.184578 | 0.1548451 | 5.13002 | 0.001 |
| Residual | 28 | 6.465508 | 0.8451549 | NA | NA |
| Total | 29 | 7.650086 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04408 0.044082 3.1744 999 0.093 .
Residuals 28 0.38883 0.013887
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.147077 | 0.1608783 | 5.368223 | 0.001 |
| Residual | 28 | 5.983012 | 0.8391217 | NA | NA |
| Total | 29 | 7.130089 | 1.0000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00007 0.0000689 0.0044 999 0.941
Residuals 28 0.43637 0.0155847 adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.1298187 | 0.1018776 | 3.17615 | 0.001 |
| Residual | 28 | 1.1444432 | 0.8981224 | NA | NA |
| Total | 29 | 1.2742619 | 1.0000000 | NA | NA |
p0<-richness %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("FMT","Transplant"),
values=c('#40714b', "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p1<-neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("FMT","Transplant"),
values=c('#40714b', "#d5b52c"))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
p2<-phylo %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(name="treatment",
breaks=c("FMT","Transplant"),
values=c('#40714b', "#d5b52c")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)
4.10.3 Structural zeros
FMT_samples <- sample_metadata %>%
filter(type == "Treatment" & treatment == "FMT") %>%
dplyr::select(sample) %>%
pull()
Transplant_samples <- sample_metadata %>%
filter(type == "Treatment" & treatment =="Transplant") %>%
dplyr::select(sample) %>% pull()
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",subset_meta$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_FMT = all(c_across(all_of(FMT_samples)) == 0)) %>%
mutate(all_zeros_Transplant = all(c_across(all_of(Transplant_samples)) == 0)) %>%
mutate(average_FMT = mean(c_across(all_of(FMT_samples)), na.rm = TRUE)) %>%
mutate(average_Transplant = mean(c_across(all_of(Transplant_samples)), na.rm = TRUE)) %>%
filter(all_zeros_FMT == TRUE || all_zeros_Transplant==TRUE) %>%
mutate(present = case_when(
all_zeros_FMT & !all_zeros_Transplant ~ "Transplant",
!all_zeros_FMT & all_zeros_Transplant ~ "FMT",
!all_zeros_FMT & !all_zeros_Transplant ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "FMT", average_FMT, average_Transplant)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Structural zeros in each group
fmt <- structural_zeros %>%
filter(present=="FMT") %>%
count(phylum, name = "FMT") %>%
arrange(desc(FMT))
fmt_f <- structural_zeros %>%
filter(present=="FMT") %>%
count(family, name = "FMT") %>%
arrange(desc(FMT)) structural_zeros %>%
filter(present=="Transplant") %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt, by="phylum" ) %>%
mutate(across(everything(), ~ replace_na(., 0))) %>%
as.data.frame() %>%
summarise(across(where(is.numeric), ~ sum(.x, na.rm = TRUE))) Transplant FMT
1 52 55Phyla to which the structural zeros belong in each group
structural_zeros %>%
filter(present=="Transplant") %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt, by="phylum" ) %>%
mutate(across(everything(), ~ replace_na(., 0))) %>%
tt()| phylum | Transplant | FMT |
|---|---|---|
| p__Bacillota_A | 20 | 21 |
| p__Bacillota | 16 | 10 |
| p__Pseudomonadota | 6 | 11 |
| p__Bacteroidota | 3 | 6 |
| p__Desulfobacterota | 2 | 2 |
| p__Bacillota_B | 1 | 0 |
| p__Chlamydiota | 1 | 0 |
| p__Cyanobacteriota | 1 | 0 |
| p__Fusobacteriota | 1 | 0 |
| p__Verrucomicrobiota | 1 | 2 |
| p__Actinomycetota | 0 | 1 |
| p__Bacillota_C | 0 | 1 |
| p__Campylobacterota | 0 | 1 |
Families to which the structural zeros belong in each group
structural_zeros %>%
filter(present=="Transplant") %>%
count(family, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt_f, by="family" ) %>%
mutate(across(everything(), ~ replace_na(., 0))) %>%
tt()| family | Transplant | FMT |
|---|---|---|
| f__Lachnospiraceae | 7 | 9 |
| f__Erysipelotrichaceae | 6 | 5 |
| f__UBA660 | 6 | 0 |
| f__Enterobacteriaceae | 5 | 2 |
| f__CAG-508 | 3 | 0 |
| f__Ruminococcaceae | 3 | 3 |
| f__Anaerovoracaceae | 2 | 0 |
| f__Coprobacillaceae | 2 | 0 |
| f__Desulfovibrionaceae | 2 | 2 |
| f__Oscillospiraceae | 2 | 1 |
| f__Tannerellaceae | 2 | 1 |
| f__UBA1242 | 2 | 0 |
| f__ | 1 | 3 |
| f__Akkermansiaceae | 1 | 0 |
| f__CAG-239 | 1 | 2 |
| f__Chlamydiaceae | 1 | 0 |
| f__Enterococcaceae | 1 | 3 |
| f__Fusobacteriaceae | 1 | 0 |
| f__Gastranaerophilaceae | 1 | 0 |
| f__Marinifilaceae | 1 | 0 |
| f__Mycoplasmoidaceae | 1 | 1 |
| f__Peptococcaceae | 1 | 0 |
| f__Anaerotignaceae | 0 | 2 |
| f__Bacteroidaceae | 0 | 2 |
| f__LL51 | 0 | 2 |
| f__UBA3700 | 0 | 2 |
| f__Acutalibacteraceae | 0 | 1 |
| f__Arcobacteraceae | 0 | 1 |
| f__CAG-274 | 0 | 1 |
| f__CALVMC01 | 0 | 1 |
| f__Devosiaceae | 0 | 1 |
| f__Mycobacteriaceae | 0 | 1 |
| f__Pumilibacteraceae | 0 | 1 |
| f__RUG11792 | 0 | 1 |
| f__Rhizobiaceae | 0 | 1 |
| f__Rikenellaceae | 0 | 1 |
| f__Sphingobacteriaceae | 0 | 1 |
| f__Streptococcaceae | 0 | 1 |
| f__UBA1997 | 0 | 1 |
| f__UBA3830 | 0 | 1 |
| f__Weeksellaceae | 0 | 1 |
4.10.3.1 Functional capacities of the structural zeros
#Aggregate bundle-level GIFTs into the compound level
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% structural_zeros$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
#Aggregate element-level GIFTs into the function level
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
#Aggregate function-level GIFTs into overall Biosynthesis, Degradation and Structural GIFTs
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
TRUE ~ treatment
))%>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
filter(genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=FMT-Transplant)%>%
mutate(group_color = ifelse(Difference <0, "Transplant","FMT")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
#scale_fill_manual() +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Elevation")
4.10.4 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("FMT1", "FMT2") ~ "FMT",
TRUE ~ treatment
))%>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant")) %>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "treatment",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output$res %>%
dplyr::select(taxon, lfc_treatmentTransplant, p_treatmentTransplant) %>%
filter(p_treatmentTransplant < 0.05) %>%
dplyr::arrange(p_treatmentTransplant) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_treatmentTransplant)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()Differential abundance MAGs in each group
fmt_count <- ancombc_rand_table_mag %>%
filter(lfc_treatmentTransplant<0) %>%
count(phylum, name = "FMT") %>%
arrange(desc(FMT))
ancombc_rand_table_mag %>%
filter(lfc_treatmentTransplant>0) %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt_count, by="phylum")%>%
mutate(across(where(is.numeric), ~ replace_na(., 0))) phylum Transplant FMT
1 Bacteroidota 13 5
2 Bacillota_A 4 17
3 Bacillota 3 1
4 Campylobacterota 1 1
5 Desulfobacterota 1 1
6 Spirochaetota 1 0
7 Pseudomonadota 0 5
8 Cyanobacteriota 0 2
9 Bacillota_B 0 1
10 Bacillota_C 0 1ancombc_rand_table_mag %>%
filter(lfc_treatmentTransplant>0) %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt_count, by="phylum")%>%
mutate(across(where(is.numeric), ~ replace_na(., 0))) %>%
as.data.frame() %>%
summarise(across(where(is.numeric), ~ sum(.x, na.rm = TRUE))) Transplant FMT
1 23 34ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_treatmentTransplant, y=forcats::fct_reorder(genome,lfc_treatmentTransplant), fill=phylum)) + #forcats::fct_rev()
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
# coord_flip()+
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))
4.10.5 Differences in the functional capacities
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
TRUE ~ treatment
))%>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.10.5.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 FMT 0.359 0.0259
2 Transplant 0.356 0.0432MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94871, p-value = 0.1561
Wilcoxon rank sum exact test
data: value by treatment
W = 128, p-value = 0.4826
alternative hypothesis: true location shift is not equal to 0element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=FMT-Transplant)%>%
mutate(group_color = ifelse(Difference <0, "Transplant","FMT")) means_gift <- element_gift_filt %>%
select(-treatment) %>%
pivot_longer(!sample, names_to = "Elements", values_to = "abundance") %>%
left_join(sample_metadata, by=join_by(sample==sample)) %>%
group_by(treatment, Elements) %>%
summarise(mean=mean(abundance))
log_fold <- means_gift %>%
group_by(Elements) %>%
summarise(
logfc_fmt_transplant = log2(mean[treatment == "FMT"] / mean[treatment == "Transplant"])
) %>%
left_join(., difference_table, by="Elements")difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=c('#40714b', "#d5b52c")) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")
log_fold %>%
filter(Elements!="D0611") %>%
ggplot(aes(x=forcats::fct_reorder(Function,logfc_fmt_transplant), y=logfc_fmt_transplant, fill=group_color)) +
geom_col() +
#scale_fill_manual() +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Log_fold")+
labs(fill="Treatment")
4.10.5.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 FMT 0.349 0.0221
2 Transplant 0.355 0.0377MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.87044, p-value = 0.001713
Wilcoxon rank sum exact test
data: value by treatment
W = 121, p-value = 0.6801
alternative hypothesis: true location shift is not equal to 0function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))
significant_functional# A tibble: 3 × 4
trait p_value p_adjust Function
<chr> <dbl> <dbl> <chr>
1 B04 0.000220 0.00441 SCFA biosynthesis
2 B10 0.00284 0.0189 Antibiotic biosynthesis
3 S02 0.00174 0.0174 Appendages